﻿using System;
using System.Diagnostics;
using System.Collections.Generic;
using System.Text;
using System.Drawing;

namespace LowLevelGraphics
{
    /// <summary>
    /// 
    /// </summary>
    [DebuggerDisplay("x = {x}, y = {y}")]
    public struct PointD
    {
        public double x;
        public double y;

        /// <summary>
        /// Initializes a new instance of the <see cref="PointD"/> struct.
        /// </summary>
        /// <param name="_x">The _x.</param>
        /// <param name="_y">The _y.</param>
        public PointD(double _x, double _y)
        {
            x = _x;
            y = _y;

        }

        /// <summary>
        /// Gets the empty.
        /// </summary>
        /// <value>The empty.</value>
        public static PointD Empty
        {
            get { return new PointD(0, 0); }
        }

        /// <summary>
        /// Ceilings the specified _point D.
        /// </summary>
        /// <param name="_pointD">The _point D.</param>
        /// <returns></returns>
        public static PointF Ceiling(PointD _pointD)
        {
            return new PointF((float)_pointD.x, (float)_pointD.y);
        }

        /// <summary>
        /// Returns a <see cref="System.String"/> that represents this instance.
        /// </summary>
        /// <returns>
        /// A <see cref="System.String"/> that represents this instance.
        /// </returns>
        public override string ToString()
        {
            return string.Format("{0,1}", x, y);
        }
    }

    /// <summary>
    /// 
    /// </summary>
    public class Maths
    {
        public static float[] gi = new float[13];
        public static float[] hi = new float[13];

        /// <summary>
        /// 
        /// </summary>
        public static float[] h = 
	    {
		    -0.00332761f,0.00569794f,	0.0196637f,	-0.0482603f,	-0.0485391f,
		    0.292562f,	0.564406f,	0.292562f,	-0.0485391f,	-0.0482602f,
		    -0.0196637f,	0.00569794f,	-0.0033276f        };

        public static float[] g = 
	    {
		    0.00332761f,  0.00569794f,	-0.0196637f,	-0.0482603f,	0.0485391f,
		    0.292562f,	-0.564406f,	0.292562f,	0.0485391f,	-0.0482602f,
		    0.0196637f,	0.00569794f,	0.0033276f        };


        /// <summary>
        /// As the specified n X.
        /// </summary>
        /// <param name="nX">The n X.</param>
        /// <param name="nXSize">Size of the n X.</param>
        /// <returns></returns>
        public static int a(int nX, int nXSize)
        {
            if (nX < 0)
                nX = -nX;
            if (nX >= nXSize)
                nX = nXSize * 2 - nX - 2;
            return nX;
        }

        /// <summary>
        /// Does a fast fourier transformation.
        /// </summary>
        /// <param name="TD">The TD.</param>
        /// <param name="FD">The FD.</param>
        /// <param name="r">The r.</param>
        public static void FastFourierTransformation(PointD[] TD, PointD[] FD, int r)
        {
            int count;
            int i, j, k, bfsize, p;
            double angle;


            /**/
            count = 1 << r;

            /**/
            //
            PointD[] W = new PointD[count / 2];
            PointD[] X1 = new PointD[count];
            PointD[] X2 = new PointD[count];
            PointD[] X = new PointD[count];

            /**/
            for (i = 0; i < count / 2; i++)
            {
                angle = -i * Math.PI * 2 / count;
                W[i].x = Math.Cos(angle);
                W[i].y = Math.Sin(angle);
            }

            /**/
            //memcpy(X1,TD,sizeof(COMPLEX)*count);
            for (i = 0; i < count; i++)
            {
                X1[i] = TD[i];
            }

            /**/
            for (k = 0; k < r; k++)
            {
                for (j = 0; j < 1 << k; j++)
                {
                    bfsize = 1 << (r - k);
                    for (i = 0; i < bfsize / 2; i++)
                    {
                        p = j * bfsize;
                        X2[i + p].x = X1[i + p].x + X1[i + p + bfsize / 2].x;
                        X2[i + p].y = X1[i + p].y + X1[i + p + bfsize / 2].y;
                        X2[i + p + bfsize / 2].x = (X1[i + p].x - X1[i + p + bfsize / 2].x) * W[i * (1 << k)].x -
                            (X1[i + p].y - X1[i + p + bfsize / 2].y) * W[i * (1 << k)].y;
                        X2[i + p + bfsize / 2].y = (X1[i + p].x - X1[i + p + bfsize / 2].x) * W[i * (1 << k)].y +
                            (X1[i + p].y - X1[i + p + bfsize / 2].y) * W[i * (1 << k)].x;
                    }
                }
                X = X1;
                X1 = X2;
                X2 = X;
            }

            /**/
            for (j = 0; j < count; j++)
            {
                p = 0;
                for (i = 0; i < r; i++)
                {
                    if ((j & (1 << i)) != 0)
                        p += 1 << (r - i - 1);
                }
                FD[j] = X1[p];
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="i"></param>
        /// <returns></returns>
        public static long sqrt(long i)
        {
            long r = 0;
            long rnew = 1;
            long rold = r;

            do
            {
                rold = r;
                r = rnew;
                rnew = (r + (i / r));
                rnew >>= 1;
            }
            while (rold != rnew);
            return rnew;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="i"></param>
        /// <returns></returns>
        public static bool IsPrime(long i)
        {
            long si, j;

            si = sqrt(i);

            for (j = 2; (j <= si); j++)
            {
                if (i % j == 0)
                    return false;
            }

            return true;
        }

        /// <summary>
        /// Fouriers the specified TD.
        /// </summary>
        /// <param name="TD">The TD.</param>
        /// <param name="lWidth">Width of the l.</param>
        /// <param name="lHeight">Height of the l.</param>
        /// <param name="FD">The FD.</param>
        /// <returns></returns>
        public static bool Fourier(PointD[] TD, int lWidth, int lHeight, PointD[] FD)
        {
            //
            long i;
            long j;
            int k;
            // 
            int w = 1;
            int h = 1;
            int wp = 0;
            int hp = 0;
            //
            while (w < lWidth)
            {
                w *= 2;
                wp++;
            }
            while (h < lHeight)
            {
                h *= 2;
                hp++;
            }
            // 
            PointD[] TempT = new PointD[h];
            PointD[] TempF = new PointD[h];
            long ji = w * h * 3;

            //
            for (i = 0; i < w * 3; i++)
            {
                // ³éÈ¡Êý¾Ý
                for (j = 0; j < h; j++)
                {
                    if (j * w * 3 + i < w * h * 3)
                        TempT[j].x = TD[j * w * 3 + i].x;
                    TempT[j].y = TD[j * w * 3 + i].y;
                }
                // Ò»Î¬¿ìËÙ¸µÁ¢Ò¶±ä»»
                FastFourierTransformation(TempT, TempF, hp);
                // 
                for (j = 0; j < h; j++)
                {
                    TD[j * w * 3 + i].x = TempF[j].x;
                    TD[j * w * 3 + i].y = TempF[j].y;
                }
            }
            //
            PointD[] TempT1 = new PointD[w];
            PointD[] TempF1 = new PointD[w];
            //
            for (i = 0; i < h; i++)
            {
                for (k = 0; k < 3; k++)
                {
                    //
                    for (j = 0; j < w; j++)
                    {
                        TempT1[j].x = TD[i * w * 3 + j * 3 + k].x;
                        TempT1[j].y = TD[i * w * 3 + j * 3 + k].y;
                    }
                    //
                    FastFourierTransformation(TempT1, TempF1, wp);
                    //
                    for (j = 0; j < w; j++)
                    {
                        FD[i * w * 3 + j * 3 + k].x = TempF1[j].x;
                        FD[i * w * 3 + j * 3 + k].y = TempF1[j].y;
                    }
                }
            }
            // ÊÍ·ÅÄÚ´æ
            return true;
        }

        /// <summary>
        /// Inverses the fast fourier transformation.
        /// </summary>
        /// <param name="FD">The FD.</param>
        /// <param name="TD">The TD.</param>
        /// <param name="power">The power.</param>
        public static void InverseFastFourierTransformation(PointD[] FD, PointD[] TD, int power)
        {
            int i, count;

            //
            count = 1 << power;

            //
            PointD[] x = new PointD[count];

            //
            //memcpy(x,FD,sizeof(COMPLEX)*count);
            for (i = 0; i < count; i++)
                x[i] = FD[i];

            //
            for (i = 0; i < count; i++)
                x[i].y = -x[i].y;

            //
            FastFourierTransformation(x, TD, power);

            /**/
            for (i = 0; i < count; i++)
            {
                TD[i].x /= count;
                TD[i].y = -TD[i].y / count;
            }
        }

        /// <summary>
        /// Is the fourier.
        /// </summary>
        /// <param name="TD">The TD.</param>
        /// <param name="lWidth">Width of the l.</param>
        /// <param name="lHeight">Height of the l.</param>
        /// <param name="FD">The FD.</param>
        public static void IFourier(PointD[] TD, int lWidth, int lHeight, PointD[] FD)
        {
            //
            long i;
            long j;
            long k;

            //
            long w = 1;
            long h = 1;
            int wp = 0;
            int hp = 0;

            //
            while (w < lWidth)
            {
                w *= 2;
                wp++;
            }

            while (h < lHeight)
            {
                h *= 2;
                hp++;
            }

            //
            //LONG lLineBytes = WIDTHBYTES(lWidth * 24);

            //

            PointD[] TempT = new PointD[w];
            PointD[] TempF = new PointD[w];

            //
            for (i = 0; i < h; i++)
            {
                for (k = 0; k < 3; k++)
                {
                    //
                    for (j = 0; j < w; j++)
                        TempF[j] = FD[i * w * 3 + j * 3 + k];

                    // 
                    InverseFastFourierTransformation(TempF, TempT, wp);

                    // ±£´æ±ä»»½á¹û
                    for (j = 0; j < w; j++)
                        FD[i * w * 3 + j * 3 + k] = TempT[j];
                }
            }

            PointD[] TempT1 = new PointD[h];
            PointD[] TempF1 = new PointD[h];

            //
            for (i = 0; i < w * 3; i++)
            {
                //
                for (j = 0; j < h; j++)
                    TempF1[j] = FD[j * w * 3 + i];

                //
                InverseFastFourierTransformation(TempF1, TempT1, hp);

                //
                for (j = 0; j < h; j++)
                    FD[j * w * 3 + i] = TempT1[j];
            }

            for (i = 0; i < h; i++)
            {
                for (j = 0; j < w * 3; j++)
                {
                    if (i < lHeight && j < 3 * w)
                        TD[i * 3 * w + j] = FD[i * w * 3 + j];
                }
            }

        }

        /// <summary>
        /// DWs the t1.
        /// </summary>
        /// <param name="lpData">The lp data.</param>
        /// <param name="nX">The n X.</param>
        /// <param name="nY">The n Y.</param>
        public static void DWT1(int[] lpData, int nX, int nY)
        {
            //
            int i, j, k, n, m;

            // 
            float temp1, temp2;
            float[] BufferX = new float[nX];
            float[] BufferY = new float[nY];

            //
            for (n = 0; n < nY; n++)
            {
                for (m = 0; m < 3; m++)
                {
                    for (i = 0; i < nX; i += 2)
                    {
                        temp1 = temp2 = 0;
                        for (j = -6; j <= 6; j++)
                            temp1 = temp1 + lpData[n * nX * 3 + Maths.a(i + j, nX) * 3 + m] * Maths.h[j + 6];
                        for (j = -5; j <= 7; j++)
                            temp2 = temp2 + lpData[n * nX * 3 + Maths.a(i + j, nX) * 3 + m] * Maths.g[j + 5];
                        BufferX[i / 2] = temp1;
                        BufferX[i / 2 + nX / 2] = temp2;
                        //MessageBox.Show(temp2.ToString());
                    }

                    // 
                    for (k = 0; k < nX; k++)
                        lpData[n * nX * 3 + k * 3 + m] = (int)BufferX[k];
                }
            }
            //
            for (n = 0; n < nX; n++)
            {
                for (m = 0; m < 3; m++)
                {
                    for (i = 0; i < nY; i += 2)
                    {
                        temp1 = temp2 = 0;
                        for (j = -6; j <= 6; j++)
                            temp1 = temp1 + lpData[Maths.a(i + j, nY) * nX * 3 + n * 3 + m] * Maths.h[j + 6];
                        for (j = -5; j <= 7; j++)
                            temp2 = temp2 + lpData[Maths.a(i + j, nY) * nX * 3 + n * 3 + m] * Maths.g[j + 5];
                        BufferY[i / 2] = temp1;
                        BufferY[i / 2 + nY / 2] = temp2;
                    }

                    // »Ø´æÊý¾Ý
                    for (k = 0; k < nY; k++)
                        lpData[k * nX * 3 + n * 3 + m] = (int)BufferY[k];
                }
            }
        }

        /*************************************************************************
        *
        * º
        *   IDWT1()
        *
        * 
        *   void	            - 
        *
        *   :
        *   
        *
        ************************************************************************/

        /// <summary>
        /// IDWs the t1.
        /// </summary>
        /// <param name="lpData">The lp data.</param>
        /// <param name="nX">The n X.</param>
        /// <param name="nY">The n Y.</param>
        public static void IDWT1(int[] lpData, int nX, int nY)
        {
            //
            int i, j, k, n, m;

            // 
            float temp1, temp2;
            float[] Buffer1 = new float[nX];
            float[] Buffer2 = new float[nY];

            for (i = 0; i < 13; i++)
            {
                Maths.hi[i] = Maths.h[12 - i];
                Maths.gi[i] = Maths.g[12 - i];
            }

            //
            for (n = 0; n < nX; n++)
            {
                for (m = 0; m < 3; m++)
                {
                    for (k = 0; k < nY / 2; k++)
                    {
                        Buffer1[k * 2] = lpData[k * nX * 3 + n * 3 + m];
                        Buffer1[k * 2 + 1] = 0;
                    }
                    for (i = 0; i < nY; i++)
                    {
                        temp1 = 0;
                        for (j = -6; j <= 6; j++)
                        {
                            temp1 = temp1 + Buffer1[Maths.a(i + j, nY)] * Maths.hi[j + 6];
                            Buffer2[i] = temp1;
                        }
                    }
                    for (k = 0; k < nY / 2; k++)
                    {
                        Buffer1[k * 2] = lpData[(k + nY / 2) * nX * 3 + n * 3 + m];
                        Buffer1[k * 2 + 1] = 0;
                    }
                    for (i = 0; i < nY; i++)
                    {
                        temp1 = 0;
                        for (j = -7; j <= 5; j++)
                            temp1 = temp1 + Buffer1[Maths.a(i + j, nY)] * Maths.gi[j + 7];
                        temp2 = temp1 + Buffer2[i];
                        Buffer2[i] = temp2;
                    }

                    // 
                    for (k = 0; k < nY; k++)
                        lpData[k * nX * 3 + n * 3 + m] = (int)Buffer2[k];
                }
            }

            //
            for (n = 0; n < nY; n++)
            {
                for (m = 0; m < 3; m++)
                {
                    for (k = 0; k < nX / 2; k++)
                    {
                        Buffer1[k * 2] = lpData[n * nX * 3 + k * 3 + m];
                        Buffer1[k * 2 + 1] = 0;
                    }

                    for (i = 0; i < nX; i++)
                    {
                        temp1 = 0;
                        for (j = -6; j <= 6; j++)
                            temp1 = temp1 + Buffer1[Maths.a(i + j, nX)] * Maths.hi[j + 6];
                        Buffer2[i] = temp1;
                    }

                    for (k = 0; k < nX / 2; k++)
                    {
                        Buffer1[k * 2] = lpData[n * nX * 3 + k * 3 + nX * 3 / 2 + m];
                        Buffer1[k * 2 + 1] = 0;
                    }

                    for (i = 0; i < nX; i++)
                    {
                        temp1 = 0;
                        for (j = -7; j <= 5; j++)
                            temp1 = temp1 + Maths.gi[j + 7] * Buffer1[Maths.a(i + j, nX)];
                        temp2 = temp1 + Buffer2[i];
                        Buffer2[i] = temp2;
                    }

                    //
                    for (k = 0; k < nX; k++)
                        lpData[n * nX * 3 + k * 3 + m] = (int)Buffer2[k] * 4;
                }
            }

        }

        /*************************************************************************
     *
     * 
     *   DCT()
     *
     * 
     *   double * f				- 
     *   double * F				- 
     *   r						
     *
     *  :
     *   
     *
     *  :
     *  
     *
     ************************************************************************/

        /// <summary>
        /// DCTs the specified f.
        /// </summary>
        /// <param name="f">The f.</param>
        /// <param name="F">The F.</param>
        /// <param name="r">The r.</param>
        public static void DCT(double[] f, double[] F, int r)
        {
            //
            int i;

            //
            double dTemp;

            //
            long N = 1 << r;

            //

            PointD[] XIn = new PointD[N * 2];
            PointD[] XOut = new PointD[N * 2];

            //
            for (i = 0; i < N * 2; i++)
            {
                XIn[i].x = 0;
                XIn[i].y = 0;
                XOut[i].x = 0;
                XOut[i].x = 0;
            }

            //
            for (i = 0; i < N; i++)
                XIn[i].x = f[i];

            //
            Maths.FastFourierTransformation(XIn, XOut, r + 1);

            //
            dTemp = 1.0f / Math.Sqrt(N);

            //
            F[0] = XOut[0].x * dTemp;

            dTemp *= Math.Sqrt(2);

            //
            for (i = 1; i < N; i++)
                F[i] = (XOut[i].x * Math.Cos(i * Math.PI / (N * 2)) + XOut[i].y * Math.Sin(i * Math.PI / (N * 2))) * dTemp;
        }
        /*************************************************************************
         *
         * º¯
         *   FreqDCT()
         *
         * 
         *   double* f			- 
         *   double* F			- 
         *	 LONG lWidth		- 
         *	 LONG lHeight		- 
         *
         * :
         *   ¸
         *
         ************************************************************************/

        /// <summary>
        /// Freqs the DCT.
        /// </summary>
        /// <param name="f">The f.</param>
        /// <param name="F">The F.</param>
        /// <param name="lWidth">Width of the l.</param>
        /// <param name="lHeight">Height of the l.</param>
        public static void FreqDCT(double[] f, double[] F, long lWidth, long lHeight)
        {
            //
            long i;
            long j;
            long k;

            //
            long w = 1;
            long h = 1;
            int wp = 0;
            int hp = 0;

            //
            while (w < lWidth)
            {
                w *= 2;
                wp++;
            }

            while (h < lHeight)
            {
                h *= 2;
                hp++;
            }

            //
            double[] TempIn = new double[h];
            double[] TempOut = new double[h];

            //
            for (i = 0; i < w * 3; i++)
            {
                //
                for (j = 0; j < h; j++)
                    TempIn[j] = f[j * w * 3 + i];

                //
                Maths.DCT(TempIn, TempOut, hp);

                //
                for (j = 0; j < h; j++)
                    f[j * w * 3 + i] = TempOut[j];
            }

            //
            double[] TempIn1 = new double[w];
            double[] TempOut1 = new double[w];

            //
            for (i = 0; i < h; i++)
            {
                for (k = 0; k < 3; k++)
                {
                    //
                    for (j = 0; j < w; j++)
                        TempIn1[j] = f[i * w * 3 + j * 3 + k];

                    //
                    DCT(TempIn1, TempOut1, wp);

                    //
                    for (j = 0; j < w; j++)
                        F[i * w * 3 + j * 3 + k] = TempOut1[j];
                }
            }
        }
        /*************************************************************************
        *
        * 
        *   IDCT()
        *
        * :
        *   double * F				- 
        *   double * f				- 
        *   r						
        *
        * :
        *   
        *
        ************************************************************************/

        public static void IDCT(double[] F, double[] f, int r)
        {
            //
            int i;

            //
            double dTemp, d0;

            //
            long N = 1 << r;

            //

            PointD[] XIn = new PointD[N * 2];
            PointD[] XOut = new PointD[N * 2];

            //
            for (i = 0; i < N * 2; i++)
            {
                XIn[i].x = 0;
                XIn[i].y = 0;
                XOut[i].x = 0;
                XOut[i].x = 0;
            }

            //
            for (i = 0; i < N; i++)
            {
                XIn[i].x = F[i] * Math.Cos(i * Math.PI / (N * 2));
                XIn[i].y = F[i] * Math.Sin(i * Math.PI / (N * 2));
            }

            //
            Maths.InverseFastFourierTransformation(XIn, XOut, r + 1);

            //
            dTemp = Math.Sqrt(2.0 / N);
            d0 = (Math.Sqrt(1.0 / N) - dTemp) * F[0];

            //
            for (i = 0; i < N; i++)
                f[i] = d0 + XOut[i].x * dTemp * 2 * N;

        }

        /// <summary>
        /// THREECROSSs the specified matrix.
        /// </summary>
        /// <param name="Matrix">The matrix.</param>
        /// <param name="QMatrix">The Q matrix.</param>
        /// <param name="MainCross">The main cross.</param>
        /// <param name="HypoCross">The hypo cross.</param>
        public static void THREECROSS(double[,] Matrix, double[,] QMatrix,
    double[] MainCross, double[] HypoCross)
        {
            //
            int i, j, k;
            double h, f, g, h2;

            //
            for (i = 0; i <= 1; i++)
                for (j = 0; j <= 1; j++)
                    QMatrix[i, j] = Matrix[i, j];

            for (i = 1; i >= 1; i--)
            {
                h = 0.0;
                if (i > 1)
                    for (k = 0; k <= i - 1; k++)
                        h = h + QMatrix[k, i] * QMatrix[k, i];

                //
                if (h + 1.0 == 1.0)
                {
                    HypoCross[i] = 0.0;
                    if (i == 1)
                        HypoCross[i] = QMatrix[i - 1, i];
                    MainCross[i] = 0.0;
                }

                    //
                else
                {
                    //
                    HypoCross[i] = Math.Sqrt(h);

                    //
                    if (QMatrix[i - 1, i] > 0.0)
                        HypoCross[i] = -HypoCross[i];

                    h = h - QMatrix[i - 1, i] * HypoCross[i];
                    QMatrix[i - 1, i] = QMatrix[i - 1, i] - HypoCross[i];
                    f = 0.0;

                    // householder
                    for (j = 0; j <= i - 1; j++)
                    {
                        QMatrix[i, j] = QMatrix[j, i] / h;
                        g = 0.0;

                        for (k = 0; k <= j; k++)
                            g = g + QMatrix[k, j] * QMatrix[k, i];

                        if (j + 1 <= i - 1)
                            for (k = j + 1; k <= i - 1; k++)
                                g = g + QMatrix[j, k] * QMatrix[k, i];

                        HypoCross[j] = g / h;
                        f = f + g * QMatrix[i, j];
                    }

                    h2 = f / (h + h);

                    //
                    for (j = 0; j <= i - 1; j++)
                    {
                        f = QMatrix[j, i];
                        g = HypoCross[j] - h2 * f;
                        HypoCross[j] = g;

                        for (k = 0; k <= j; k++)
                            QMatrix[k, j] = QMatrix[k, j] - f * HypoCross[k] - g * QMatrix[k, i];
                    }
                    MainCross[i] = h;
                }
            }

            // ¸³ÁãÖµ
            for (i = 0; i <= 0; i++)
                HypoCross[i] = HypoCross[i + 1];
            HypoCross[1] = 0.0;
            MainCross[0] = 0.0;

            for (i = 0; i <= 1; i++)
            {
                //
                if ((MainCross[i] != 0.0) && (i - 1 >= 0))
                    for (j = 0; j <= i - 1; j++)
                    {
                        g = 0.0;
                        for (k = 0; k <= i - 1; k++)
                            g = g + QMatrix[k, i] * QMatrix[j, k];
                        for (k = 0; k <= i - 1; k++)
                            QMatrix[j, k] = QMatrix[j, k] - g * QMatrix[i, k];
                    }

                //
                MainCross[i] = QMatrix[i, i];
                QMatrix[i, i] = 1.0;

                //
                if (i - 1 >= 0)
                    for (j = 0; j <= i - 1; j++)
                    {
                        QMatrix[j, i] = 0.0;
                        QMatrix[i, j] = 0.0;
                    }
            }

        }

        /// <summary>
        /// Discretes the wavelet transform.
        /// </summary>
        /// <param name="lpData">The lp data.</param>
        /// <param name="nX">The n X.</param>
        /// <param name="nY">The n Y.</param>
        public static void DiscreteWaveletTransform(int[] lpData, int nX, int nY)
        {
            // 
            int i, j, k, n;

            // 
            float temp1, temp2;
            float[] BufferX = new float[nX];
            float[] BufferY = new float[nY];

            //
            for (n = 0; n < nY; n++)
            {
                for (i = 0; i < nX; i += 2)
                {
                    temp1 = temp2 = 0;
                    for (j = -6; j <= 6; j++)
                        temp1 = temp1 + lpData[n * nX + a(i + j, nX)] * h[j + 6];
                    for (j = -5; j <= 7; j++)
                        temp2 = temp2 + lpData[n * nX + a(i + j, nX)] * g[j + 5];
                    BufferX[i / 2] = temp1;
                    BufferX[i / 2 + nX / 2] = temp2;
                    //MessageBox.Show(temp2.ToString());
                }

                //
                for (k = 0; k < nX; k++)
                    lpData[n * nX + k] = (int)BufferX[k];
            }

            //
            for (n = 0; n < nX; n++)
            {
                for (i = 0; i < nY; i += 2)
                {
                    temp1 = temp2 = 0;
                    for (j = -6; j <= 6; j++)
                        temp1 = temp1 + lpData[a(i + j, nY) * nX + n] * h[j + 6];
                    for (j = -5; j <= 7; j++)
                        temp2 = temp2 + lpData[a(i + j, nY) * nX + n] * g[j + 5];
                    BufferY[i / 2] = temp1;
                    BufferY[i / 2 + nY / 2] = temp2;
                }

                // »Ø´æÊý¾Ý
                for (k = 0; k < nY; k++)
                    lpData[k * nX + n] = (int)BufferY[k];
            }
        }

        /// <summary>
        /// Inverses the dicrete wavelet transform.
        /// </summary>
        /// <param name="lpData">The lp data.</param>
        /// <param name="nX">The n X.</param>
        /// <param name="nY">The n Y.</param>
        public static void InverseDicreteWaveletTransform(int[] lpData, int nX, int nY)
        {
            // Ñ­»·±äÁ¿
            int i, j, k, n;

            // ÖÐ¼ä±äÁ¿
            float temp1, temp2;
            float[] Buffer1 = new float[nX];
            float[] Buffer2 = new float[nY];

            for (i = 0; i < 13; i++)
            {
                hi[i] = h[12 - i];
                gi[i] = g[12 - i];
            }

            // ´¹Ö±·½Ïò
            for (n = 0; n < nX; n++)
            {
                for (k = 0; k < nY / 2; k++)
                {
                    Buffer1[k * 2] = lpData[k * nX + n];
                    Buffer1[k * 2 + 1] = 0;
                }
                for (i = 0; i < nY; i++)
                {
                    temp1 = 0;
                    for (j = -6; j <= 6; j++)
                    {
                        temp1 = temp1 + Buffer1[a(i + j, nY)] * hi[j + 6];
                        Buffer2[i] = temp1;
                    }
                }
                for (k = 0; k < nY / 2; k++)
                {
                    Buffer1[k * 2] = lpData[(k + nY / 2) * nX + n];
                    Buffer1[k * 2 + 1] = 0;
                }
                for (i = 0; i < nY; i++)
                {
                    temp1 = 0;
                    for (j = -7; j <= 5; j++)
                        temp1 = temp1 + Buffer1[a(i + j, nY)] * gi[j + 7];
                    temp2 = temp1 + Buffer2[i];
                    Buffer2[i] = temp2;
                }

                //
                for (k = 0; k < nY; k++)
                    lpData[k * nX + n] = (int)Buffer2[k];
            }

            //
            for (n = 0; n < nY; n++)
            {
                for (k = 0; k < nX / 2; k++)
                {
                    Buffer1[k * 2] = lpData[n * nX + k];
                    Buffer1[k * 2 + 1] = 0;
                }

                for (i = 0; i < nX; i++)
                {
                    temp1 = 0;
                    for (j = -6; j <= 6; j++)
                        temp1 = temp1 + Buffer1[a(i + j, nX)] * hi[j + 6];
                    Buffer2[i] = temp1;
                }

                for (k = 0; k < nX / 2; k++)
                {
                    Buffer1[k * 2] = lpData[n * nX + k + nX / 2];
                    Buffer1[k * 2 + 1] = 0;
                }

                for (i = 0; i < nX; i++)
                {
                    temp1 = 0;
                    for (j = -7; j <= 5; j++)
                        temp1 = temp1 + gi[j + 7] * Buffer1[a(i + j, nX)];
                    temp2 = temp1 + Buffer2[i];
                    Buffer2[i] = temp2;
                }

                //
                for (k = 0; k < nX; k++)
                    lpData[n * nX + k] = (int)Buffer2[k] * 4;
            }

        }

        /// <summary>
        /// BSTQs the specified main cross.
        /// </summary>
        /// <param name="MainCross">The main cross.</param>
        /// <param name="HypoCross">The hypo cross.</param>
        /// <param name="Matrix">The matrix.</param>
        /// <param name="Eps">The eps.</param>
        /// <param name="MaxT">The max T.</param>
        public static void BSTQ(double[] MainCross, double[] HypoCross,
            double[,] Matrix, double Eps, int MaxT)
        {
            //
            int i, j, k, m, it;
            double d, f, h, g, p, r, e, s;

            //
            HypoCross[1] = 0.0;
            d = 0.0;
            f = 0.0;

            for (j = 0; j <= 1; j++)
            {
                // 
                it = 0;
                h = Eps * (Math.Abs(MainCross[j]) + Math.Abs(HypoCross[j]));
                if (h > d)
                    d = h;
                m = j;

                while ((m <= 1) && (Math.Abs(HypoCross[m]) > d))
                    m = m + 1;

                if (m != j)
                {
                    //
                    do
                    {
                        //
                        if (it == MaxT)
                            return;
                        it = it + 1;
                        g = MainCross[j];
                        p = (MainCross[j + 1] - g) / (2.0 * HypoCross[j]);
                        r = Math.Sqrt(p * p + 1.0);

                        //
                        if (p >= 0.0)
                            MainCross[j] = HypoCross[j] / (p + r);
                        else
                            MainCross[j] = HypoCross[j] / (p - r);

                        h = g - MainCross[j];

                        //
                        for (i = j + 1; i <= 1; i++)
                            MainCross[i] = MainCross[i] - h;

                        //
                        f = f + h;
                        p = MainCross[m];
                        e = 1.0; s = 0.0;

                        for (i = m - 1; i >= j; i--)
                        {
                            g = e * HypoCross[i];
                            h = e * p;

                            //
                            if (Math.Abs(p) >= Math.Abs(HypoCross[i]))
                            {
                                e = HypoCross[i] / p;
                                r = Math.Sqrt(e * e + 1.0);
                                HypoCross[i + 1] = s * p * r;
                                s = e / r; e = 1.0 / r;
                            }
                            else
                            {
                                e = p / HypoCross[i];
                                r = Math.Sqrt(e * e + 1.0);
                                HypoCross[i + 1] = s * HypoCross[i] * r;
                                s = 1.0 / r; e = e / r;
                            }

                            p = e * MainCross[i] - s * g;
                            MainCross[i + 1] = h + s * (e * g + s * MainCross[i]);

                            //
                            for (k = 0; k <= 1; k++)
                            {
                                h = Matrix[i + 1, k];
                                Matrix[i + 1, k] = s * Matrix[i, k] + e * h;
                                Matrix[i, k] = e * Matrix[i, k] - s * h;
                            }

                        }

                        //
                        HypoCross[j] = s * p;
                        MainCross[j] = e * p;

                    }
                    while (Math.Abs(HypoCross[j]) > d);
                }

                MainCross[j] = MainCross[j] + f;
            }

            //
            for (i = 0; i <= 1; i++)
            {
                k = i; p = MainCross[i];

                //
                if (i + 1 <= 2 - 1)
                {
                    j = i + 1;
                    while ((j <= 1) && (MainCross[j] <= p))
                    {
                        k = j;
                        p = MainCross[j];
                        j = j + 1;
                    }
                }

                //
                if (k != i)
                {
                    MainCross[k] = MainCross[i];
                    MainCross[i] = p;
                    for (j = 0; j <= 1; j++)
                    {
                        p = Matrix[i, j];
                        Matrix[i, j] = Matrix[k, j];
                        Matrix[k, j] = p;
                    }
                }
            }
        }

        /// <summary>
        /// Is the freq DCT.
        /// </summary>
        /// <param name="f">The f.</param>
        /// <param name="F">The F.</param>
        /// <param name="lWidth">Width of the l.</param>
        /// <param name="lHeight">Height of the l.</param>
        public static void IFreqDCT(double[] f, double[] F, long lWidth, long lHeight)
        {
            //
            long i;
            long j;
            long k;

            //
            long w = 1;
            long h = 1;
            int wp = 0;
            int hp = 0;

            //
            while (w < lWidth)
            {
                w *= 2;
                wp++;
            }

            while (h < lHeight)
            {
                h *= 2;
                hp++;
            }

            //
            long lLineBytes = lWidth * 3;

            //
            double[] TempIn = new double[w];
            double[] TempOut = new double[w];

            //
            for (i = 0; i < h; i++)
            {
                for (k = 0; k < 3; k++)
                {
                    //
                    for (j = 0; j < w; j++)
                        TempIn[j] = F[i * w * 3 + j * 3 + k];

                    //
                    IDCT(TempIn, TempOut, wp);

                    //
                    for (j = 0; j < w; j++)
                        F[i * w * 3 + j * 3 + k] = TempOut[j];
                }
            }

            double[] TempIn1 = new double[h];
            double[] TempOut1 = new double[h];

            //
            for (i = 0; i < w * 3; i++)
            {
                // 
                for (j = 0; j < h; j++)
                    TempIn1[j] = F[j * w * 3 + i];

                //
                IDCT(TempIn1, TempOut1, hp);

                //
                for (j = 0; j < h; j++)
                    F[j * w * 3 + i] = TempOut1[j];
            }


            for (i = 0; i < h; i++)
            {
                for (j = 0; j < w * 3; j++)
                {
                    if (i < lHeight && j < lLineBytes)
                        f[i * lLineBytes + j] = F[i * w * 3 + j];
                }
            }
        }

        /// <summary>
        /// Freqs the WALSH.
        /// </summary>
        /// <param name="f">The f.</param>
        /// <param name="F">The F.</param>
        /// <param name="lWidth">Width of the l.</param>
        /// <param name="lHeight">Height of the l.</param>
        public static void FreqWALSH(double[] f, double[] F, int lWidth, int lHeight)
        {
            //
            int i;
            int j;
            int k;

            //
            int w = 1;
            int h = 1;
            int wp = 0;
            int hp = 0;

            //
            while (w < lWidth)
            {
                w *= 2;
                wp++;
            }

            while (h < lHeight)
            {
                h *= 2;
                hp++;
            }

            //
            double[] TempIn = new double[h];
            double[] TempOut = new double[h];

            //
            //
            for (i = 0; i < w * 3; i++)
            {
                //
                for (j = 0; j < h; j++)
                    TempIn[j] = f[j * w * 3 + i];

                //
                WALSH(TempIn, TempOut, hp);


                for (j = 0; j < h; j++)
                    f[j * w * 3 + i] = TempOut[j];
            }

            //
            double[] TempIn1 = new double[w];
            double[] TempOut1 = new double[w];

            //
            for (i = 0; i < h; i++)
            {
                for (k = 0; k < 3; k++)
                {
                    //
                    for (j = 0; j < w; j++)
                        TempIn1[j] = f[i * w * 3 + j * 3 + k];

                    // 
                    WALSH(TempIn1, TempOut1, wp);

                    //
                    for (j = 0; j < w; j++)
                        F[i * w * 3 + j * 3 + k] = TempOut1[j];
                }
            }
        }

        /// <summary>
        /// WALSHs the specified f.
        /// </summary>
        /// <param name="f">The f.</param>
        /// <param name="F">The F.</param>
        /// <param name="r">The r.</param>
        public static void WALSH(double[] f, double[] F, int r)
        {
            //
            int i;
            int j;
            int k;

            //
            int p;

            //
            long N = 1 << r;

            //
            double[] X = new double[N];
            double[] X1 = new double[N];
            double[] X2 = new double[N];

            // 
            //memcpy(X1, f, sizeof(double) * N);
            for (i = 0; i < N; i++)
                X1[i] = f[i];

            // 
            for (k = 0; k < r; k++)
            {
                for (j = 0; j < 1 << k; j++)
                {
                    for (i = 0; i < 1 << (r - k - 1); i++)
                    {
                        p = j * (1 << (r - k));
                        X2[i + p] = X1[i + p] + X1[i + p + (int)(1 << (r - k - 1))];
                        X2[i + p + (int)(1 << (r - k - 1))] = X1[i + p] - X1[i + p + (int)(1 << (r - k - 1))];
                    }
                }

                // 
                X = X1;
                X1 = X2;
                X2 = X;
            }

            // 
            for (j = 0; j < N; j++)
            {
                p = 0;
                for (i = 0; i < r; i++)
                {
                    if ((j & (1 << i)) != 0)
                    {
                        p += 1 << (r - i - 1);
                    }
                }

                F[j] = X1[p] / N;
            }
        }

        /// <summary>
        /// IWALSHs the specified F.
        /// </summary>
        /// <param name="F">The F.</param>
        /// <param name="f">The f.</param>
        /// <param name="r">The r.</param>
        public static void IWALSH(double[] F, double[] f, int r)
        {
            //
            int i;

            //
            long N = 1 << r;

            //
            WALSH(F, f, r);

            //
            for (i = 0; i < N; i++)
                f[i] *= N;
        }

        /// <summary>
        /// Is the freq WALSH.
        /// </summary>
        /// <param name="f">The f.</param>
        /// <param name="F">The F.</param>
        /// <param name="lWidth">Width of the l.</param>
        /// <param name="lHeight">Height of the l.</param>
        public static void IFreqWALSH(double[] f, double[] F, long lWidth, long lHeight)
        {
            //
            long i;
            long j;
            long k;

            // 
            long w = 1;
            long h = 1;
            int wp = 0;
            int hp = 0;

            //
            while (w < lWidth)
            {
                w *= 2;
                wp++;
            }

            while (h < lHeight)
            {
                h *= 2;
                hp++;
            }

            //
            long lLineBytes = lWidth * 3;

            // 
            double[] TempIn = new double[w];
            double[] TempOut = new double[w];

            //
            for (i = 0; i < h; i++)
            {
                for (k = 0; k < 3; k++)
                {
                    // 
                    for (j = 0; j < w; j++)
                        TempIn[j] = F[i * w * 3 + j * 3 + k];


                    IWALSH(TempIn, TempOut, wp);

                    //
                    for (j = 0; j < w; j++)
                        F[i * w * 3 + j * 3 + k] = TempOut[j];
                }
            }

            double[] TempIn1 = new double[h];
            double[] TempOut1 = new double[h];

            // 
            for (i = 0; i < w * 3; i++)
            {
                // 
                for (j = 0; j < h; j++)
                    TempIn1[j] = F[j * w * 3 + i];


                IWALSH(TempIn1, TempOut1, hp);

                //
                for (j = 0; j < h; j++)
                    F[j * w * 3 + i] = TempOut1[j];
            }


            for (i = 0; i < h; i++)
            {
                for (j = 0; j < w * 3; j++)
                {
                    if (i < lHeight && j < lLineBytes)
                        f[i * lLineBytes + j] = F[i * w * 3 + j];
                }
            }
        }
    }
}
